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ABSTRACT 


A  comprehensive  equation  of  state  based  on  the 
Helmholtz  free  energy  has  been  developed  for  representing 
the  thermodynamic  properties  of  pure  substances  over  a 
wide  range  of  conditions  including  the  region  near  the 
liquid-vapor  critical  point.  The  formulation  comprises 
a  set  of  terms  that  are  regular  (expandable  in  a  Taylor 
series)  in  the  vicinity  of  this  point  plus  a  term  which 
becomes  singular  there.  The  latter  takes  accoxont  of 
the  anomalous  behavior  observed  when  the  critical  state 
is  approached,  and  is  needed  to  attain  a  consistent, 
quantitative  approximation  to  experimental  data. 

Both  the  singular  and  the  nonsingular  parts  of  the 
free  energy  are  given  in  closed  form  as  functions  of 
density  and  temperature.  Moreover,  the  expressions  are 
designed  to  vary  in  an  appropriate  fashion  when  one  or 
the  other  of  the  arguments  approaches  zero  or  infinity. 
Thus,  the  limiting  form  is  in  accord  with  Debye ^s  T^  law 
for  specific  heats  at  low  temperatures,  while  it  reduces 
to  the  ideal  gas  law  when  the  density  is  allowed  to  vanish. 
The  majority  of  the  parameters  occur  as  coefficients  of 
the  various  terms,  which  facilitates  a  least-squares 
determination  of  their  values.  The  few  that  are 
involved  in  a  nonlinear  manner  have  been  introduced  to 
characterize  behavior  peculiar  to  certain  regions. 

The  work  described  here  augments  that  previously 
documented  in  Naval  Ship  Research  and  Development  Center 
Report  3664  (May  1971) .  A  more  explicit  treatment  of  the 
critical  region  is  presented  and  the  approximation  for 
the  singular  contribution  is  extended  to  a  higher  order 
to  permit  better  agreement  with  observation. 


ADMINISTRATIVE  INFORMATION 

The  work  reported  here  was  carried  out  within  the  Advanced  Software 
R&D  Group  of  the  Computer  Sciences  Division  under  the  sponsorship  of 
NAV SHIPS,  Task  Area  SF53532107,  Work  Unit  1-1830-005. 

This  material  was  originally  presented  at  the  74th  National  Meeting 
of  the  American  Institute  of  Chemical  Engineers  held  in  New  Orleans  in 
March  1973.  The  format  follows  the  guidelines  prescribed  for  papers 
presented  at  that  meeting. 
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SUMMARY 


A  fundamental  equation  of  state  based  on  the  Helmholtz  free 
energy  is  given  which  permits  correlation  of  equilibrium 
thermodynamic  properties  of  simple  fluids  over  a  wide  range 
of  conditions  embracing  both  the  dispersed  and  condensed 
phases.  The  formulation  takes  account  of  the  singularity 
at  the  critical  point  while  exhibiting  the  proper  asymptotic 
behavior  for  extreme  values  of  density  and  temperature. 

I.  INTRODUCTION 

This  paper  deals  with  the  mathematical  representation  of  the 
thermodynamic  properties  of  a  chemically  homogeneous,  stable  molecular 
substance.  It  is  intended  for  applications  in  which  a  single,  compact 
formulation  yielding  accurate  values  over  a  wide  range  of  conditions  is 
required. 

There  is  a  growing  awareness  that  comprehensive,  computer-oriented 
descriptions  of  material  properties  can  be  of  great  value  to  such  tasks 
as  engineering  design  and  contract  preparation,  as  well  as  to  the  quest 
for  new  knowledge.  Such  needs  were  recognized  long  ago  in  the  case  of 
water  due  to  its  nearly  universal  employment  as  the  working  fluid  in 
power-generating  machinery.  Since  1921,  under  the  sponsorship  of  what 
is  now  called  the  International  Association  for  the  Properties  of  Steam, 
there  has  been  a  continuing  effort  to  reassess  and  improve  the  accuracy 
of  available  data  and  to  extend  the  range  of  observation  so  that  reliable 
tables  and  formulas  can  be  compiled  for  designing  boilers  and  turbines. 

The  incentive  for  undertaking  the  study  reported  here  grew  out  of 
an  attempt  to  devise  an  effective  interpolation  scheme  for  the  data  com¬ 
prising  the  1963  International  Skeleton  Table  for  Compressed  Water  and 
Superheated  Steam.  This  set  of  values  was  adopted  by  the  Sixth  International 
Conference  on  the  Properties  of  Steam  (ref.  1) ,  and  a  concerted  effort  was 
then  made  by  groups  in  several  countries  to  produce  a  suitable  formulation. 
However,  because  of  the  great  range  covered  by  the  table  (from  1  to  1000 
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atmospheres  and  0  to  800  degrees  C) ,  no  single  equation  could  be  found 
that  would  satisfy  the  stringent  tolerances  prescribed  by  the  Sixth 
ICPS  .  An  interim  solution  was  eventually  agreed  upon  wherein  the  data 
points  were  divided  into  six  subregions,  each  having  its  own  formulation, 
with  the  stipulation  that  they  should  closely  conform  at  their  common 
boundaries . 

The  representation  finally  adopted  contained  a  total  of  151  constants 
plus  a  large  assortment  of  integer  as  well  as  fractional  exponents  (ref,  4), 
Such  a  solution  leaves  much  to  be  desired,  both  practically  and  aesthetically. 
It  became  apparent  that  a  more  rational  approach  would  have  to  be  pursued 
to  achieve  any  appreciable  simplification  of  the  situation. 

II.  OBJECTIVES 

The  notion  of  the  continuity  of  state  for  a  given  substance  implies  the 
existence  of  a  single  analytic  function  that  will  completely  characterize 
its  equilibrium  thermodynamic  properties  (ref.  5).  The  problem  then  reduces 
to  that  of  devising  an  accurate  and  concise  approximation  for  this  function. 

In  view  of  the  remarkably  similar  behavior  of  nearly  all  substances,  it 
seemed  desirable  that  such  a  formulation  should  expressly  reflect  those 
attributes  held  in  common,  so  that  it  could  have  general  application  rather 
than  be  tailored  to  a  particular  substance.  However,  it  was  decided 
to  restrict  the  investigation  to  a  description  of  the  behavior  of  pure 
substances,  and  thereby  avoid  the  added  complications  that  arise  in  dealing 
with  mixtures.  Logic  dictates  that  an  acceptable  means  of  treating  simple 
fluids  should  be  worked  out  first. 


Two  of  the  papers  presented  at  the  Seventh  ICPS  in  1968  came  close  to 
satisfying  the  requirements  (viz. ,  reproduction  of  density  values  to 
better  than  0.1  percent  over  most  of  the  region),  one  by  Keenan,  Keyes, 
Hill  and  Moore  (ref.  3)  and  the  other  by  the  present  author  (later  ex¬ 
panded  into  ref.  2).  Both  these  formulations  were  based  on  the  Helmholtz 
function,  but  each  suffered  from  the  drawback  that  no  provision  had  been 
made  for  the  singular  behavior  of  this  function  at  the  critical  point. 
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A  number  of  considerations,  both  practical  and  theoretical,  entered 
into  the  step-by-step  development  of  a  comprehensive  equation  of  state. 
Experimental  measurements,  by  themselves,  do  not  usually  fulfill  all  the 
requirements  of  users,  since  observation  is  limited  by 

1)  the  availability  of  suitable  instrumentation 

2)  the  range  over  which  measurement  is  feasible 

3)  the  precision  attainable  in  each  region 

4)  the  difficulty  and  cost  of  experimentation. 

A  mathematical  formulation  helps  to  mitigate  some  of  these  restrictions 
while  providing  for  a  variety  of  user  needs,  such  as 

1)  smoothing  and  interpolation  of  data  values 

2)  correlation  of  diverse  properties 

3)  compact  representation  of  extensive  observations 

4)  evaluation  of  quantities  not  directly  measurable,  e.g.,  entropy. 

Lastly,  there  are  constraints  which  derive  from  the  physical  behavior  of 
material  substances.  Among  the  phenomena  that  must  be  taken  into  account 
are 

1)  the  nearly  ideal  behavior  at  low  densities 

2)  the  vapor  pressure  at  low  temperatures 

3)  the  Joule-Thomson  coefficient 

4)  the  specific  heat  of  the  condensed  phases 

5)  the  specific-heat  anomaly  at  the  critical  point 

6)  the  vapor-liquid  equilibrium  curve. 

In  addition,  at  every  point  for  which  the  substance  exhibits  either  stable 
or  metastable  behavior,  a  proper  formulation  must  satisfy  the  Thermodynamic 
Inequalities,  i.e.,  both  the  isochoric  specific  heat  and  the  isothermal 
compressibility  must  be  positive.  All  in  all,  a  general  phenomenological 
description  consistent  with  established  theory  appeared  to  be  a  reasonable 
goal. 
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Ill .  OVERALL  APPROACH 


It  was  concluded  that  the  various  objectives  could  best  be 
achieved  through  the  adoption  of  one  of  the  canonical  thermodynamic 
functions.  All  equilibrium  properties  can  then  be  expressed  in  terms  of 
the  set  of  its  partial  derivatives  with  respect  to  two  appropriate 
conjugate  variables.  The  Helmholtz  free  energy  was  chosen  not  only 
because  of  its  basic  role  in  statistical  physics  but  also  because  of  its 
convenient  attributes  and  regular  behavior.  It  is  the  only  one  of  the 
canonical  functions  which  involves  directly  measurable  quantities  (such 
as  pressure,  volume,  temperature,  and  energy  in  the  form  of  heat  or  work) 
and  which  also  has  the  property  that  all  its  first-  and  second-order 
differential  coefficients  are  finite  and  continuous  for  all  physically 
realizable  values  of  the  two  independent  variables,  excluding  only  the 
critical  point. 

This  function  may  be  conveniently  used  in  either  of  two  forms: 


f  =  f(v,T)  [free  energy  per  xmit  mass] 

a  ^  a(p,T)  [free  energy  per  unit  volume] 


where  v,  p,  and  T  stand  for  the  specific  volume,  density,  and  absolute 
temperature,  respectively.  The  relation  between  the  two  forms  is  simply 
a  =  pf.  Quantities  such  as  pressure  p,  isothermal  compressibility  k^, 
specific  entropy  s,  and  specific  heat  at  constant  volume  c^  are  directly 
derivable  from  f(v,T).  Thus  we  have 


p  =  -(9f/9v) 

1/k^  =  v(92f/9v2) 

s  =  -(9f/9T) 

=  -T(92f/3T2) 


(1) 

(2) 

(3) 

(4) 


i 

I 
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other  properties  are  readily  obtained  from  various  combinations  of 
the  foregoing,  e.g., 


u 

=  f 

+ 

Ts 

(5) 

g 

=  f 

+ 

pv 

(6) 

h 

=  u 

+ 

pv  =  g  +  Ts 

(7) 

in  which  u  denotes  the  internal  energy  per  unit  mass,  g,  the  Gibbs 
function  per  unit  mass,  and  h,  the  enthalpy  per  unit  mass. 

One  of  the  main  conclusions  reached  in  the  course  of  this  study 
was  that  the  free  energy  must  be  resolved  into  a  singular  part  which 
predominates  in  the  critical  region  plus  a  regular  part  which  suffices 
elsewhere.  This  conclusion  is  based  on  experimental  investigations  of 
many  different  gases  carried  out  by  Voronel’  and  his  co-workers  in  the 
Soviet  Union  during  the  1960’s.  Their  measurements  reveal  that  the 
isochoric  specific  heat  of  pure  substances  becomes  indefinitely  large  as 
the  critical  state  is  approached  (ref.  6).  No  function  which  is  analytic 
(i.e.,  expandable  in  a  Taylor  series)  at  the  critical  point  can  adequately 
represent  behavior  of  this  type.  However,  such  a  function  can  serve  as 
the  background  onto  which  a  properly  formulated  anomaly  may  be  super¬ 
imposed.  We  shall  consider  the  regular  part  first,  and  then  introduce 
an  expression  for  the  singular  contribution. 

IV.  REGULAR  TERMS  OF  THE  APPROXIMATION 

It  will  be  convenient  in  what  follows  to  employ  nondimensional 
independent  variables,  namely. 


p/Pc 

[reduced  density] 

(8) 

T/T 

c 

[reduced  temperature] 

(9) 

where  and  T^  stand  for  the  critical  density  and  temperature  of  the 
substance  in  question.  We  proceed  by  defining  an  auxiliary  function  of 
these  two  variables  which  exhibits  a  marked  change  of  character  in  the 
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transition  from  a  dispersed  to  a  condensed  state.  Accordingly,  we  write 


<P(^,e) 


as  g  0  or  9  “ 


(10) 


=  1  asg->-“or0-)“O 

and  in  which  ri(9)  =  tio(1  “  e  (H) 

cCe)  =  eV  (12)  ‘ 

O 

n  is  the  number  of  rotational  degrees  of  freedom 
of  the  molecule 


Q  stands  for  the  molecular  partition  function  excluding 
the  factors  for  translation  and  rotation 

A 

and  e,  no  9  adjustable  constants.  A  significant  property  of  ())(5,0) 
is  that,  although  its  value  is  quite  small  under  conditions  approaching 
those  of  an  ideal  gas,  it  tends  toward  unity  at  low  temperatures  and/or 
very  high  densities.  This  feature  is  crucial  to  giving  the  subsequent 
formulation  a  wide  scope.  The  function  n(0)  is  introduced  to  render  a 
better  description  of  the  Joule-Thomson  coefficient  as  well  as  the  second 
vlrial  coefficient  at  low  temperatures.  On  the  other  hand,  C(0)  relates 
strictly  to  the  behavior  of  the  ideal  gas  and  is  assumed,  in  the  context 
of  this  work,  to  be  a  known  function.  Moreover,  by  explicitly  isolating 
the  factor  0^  from  ?(0),  we  ensure  conformity  with  Debye^s  formula  for  the 
specific  heat  of  the  condensed  phases  at  low  temperatures. 

The  different  classes  of  terms  that  contribute  to  the  regular  part 
of  the  equation  for  the  free  energy  are 

a)  the  ideal  gas  approximation 

|rT  In  ())  (13) 

b)  the  residual  entropy  at  absolute  zero 

R  In  r2o 


(14) 


(15) 


c)  the  internal  energy  at  absolute  zero 

^  p 

RT  y  A,  r 

d)  the  temperature-dependent,  hlgh-denslty  corrections 

m  I 


(16) 


in  which  R  is  the  gas  constant,  Qq  stands  for  the  (effective)  number  of 
distinct  orientations  that  the  molecules  can  assume  in  the  formation  of 
the  crystalline  phase  of  the  substance,  and  are  adjustable  coef¬ 
ficients,  and  m  is  an  index  which  fixes  the  degree  of  the  polynomials 
(15)  and  (16).  This  index  is  set  to  provide  as  many  terms  as  are  necessary 
to  fit  the  data  to  within  the  prescribed  tolerances.  It  should  be  noted 
that  none  of  the  properties  of  ideal  gases  are  in  any  way  affected  by 
inclusion  of  the  various  supplemental  terms. 

A  discussion  of  the  rationale  behind  the  design  of  the  above 
expressions,  together  with  a  detailed  analysis  of  their  as3niiptotic 
behavior,  is  given  in  reference  2.  It  will  suffice  here  to  indicate 
the  limiting  forms  to  which  the  proposed  equation  of  state  reduces  in 
each  of  four  extreme  cases : 

1)  Debye’s  law  for  solids  (as  T  0  for  constant  p) 

2)  van  der  Waals’  equation  (as  T  o®  for  constant  p) 

3)  the  ideal-gas  law  (as  p  0  for  constant  T) 

4)  a  polynomial  in  density  (as  p  ->  <»  for  constant  T). 

A  check  on  the  adequacy  of  such  an  equation  was  carried  out  some 
time  ago  on  the  basic  set  of  values  for  water  and  steam  given  in  the 
1963  International  Skeleton  Tables  over  a  grid  comprising  580  combinations 
of  pressure  and  temperature.  Both  enthalpy  and  PVT  data  were  used  as 
input  to  the  fitting  procedure.  The  resulting  PVT  approximation  (with 
the  index  m  equal  to  5)  remained  generally  within  the  assigned  tolerances 
and  was  about  as  accurate  as  that  obtained  by  Keenan  et  al  (ref.  3). 
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However,  the  concomitant  approximation  for  enthalpy  did  not  prove  to 
be  as  satisfactory,  particularly  in  the  area  close  to  the  critical 
point.  The  deficiencies  were  most  evident  in  values  calculated  for  the 
isochoric  specific  heat,  which  deviated  quite  significantly  from 
experimental  data  taken  in  this  region  (ref.  7).  This  observation 
served  as  the  motivation  for  the  next  stage  of  development. 

V.  THE  CRITICAL  REGION 

Attempts  to  extend  analytic  formulations  to  include  the  area 
around  the  critical  point  have  generally  encountered  major  obstacles. 
Either  an  inordinate  number  of  terms  must  be  used,  having — for  the 
most  part — rather  large  numerical  coefficients  but  nearly  canceling 
one  another  (which  can  lead  to  erratic  derivatives) ,  or  else  the 
tolerances  have  to  be  greatly  relaxed.  An  accuracy  of  0.1  percent 
or  less  in  fitting  pressure  (which  is  one  of  the  least  sensitive 
variables  in  the  critical  region)  cannot  be  achieved  without  some 
recognition  of  the  peculiar  properties  of  a  fluid  near  its  critical 
state  (see  Figures  1  and  2).  Careful  observation  has  revealed  striking 
departures  from  the  classical  behavior  exemplified  by  an  equation  such 
as  that  of  van  der  Waals  (ref.  8). 

The  significant  phenomena  are:  a)  the  greatly  diminished  curvature 
of  the  coexistence  or  vapory-liquid  equilibrium  curve  near  the  critical 
point;  b)  the  pronounced  flatness  of  the  pressure  isotherms  in  the 
region  lying  just  above  this  point;  and  c)  the  presence  of  an  anomaly 
in  the  specific  heat  which,  in  the  same  region,  varies  approximately  as 

-c^  ^  In  T  =  ln(T/T^-l)  for  p=^p^,  T>T^.  (17) 

It  has  been  noted  by  Yang  and  Yang  (ref.  9)  that  the  thermodynamic 
identity  relating  the  specific  heat,  pressure,  and  chemical  potential 
(or  Gibbs  function) 

pc^  =  T(92p/3T2)^  -  pT02g/3T2)^  (18) 
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Figure  1 


Shape  of  Isotherms  in  Critical  Regior 


UJ  cc  < 

c/5  CD  Q 


implies  that  one  (or  both)  of  these  second  partial  derivatives  entails 
the  same  singularity  as  the  specific  heat. 

The  hypothesis  that  the  chemical  potential  is  analytic  at  the 
critical  point  has  been  advanced  by  a  number  of  investigators  on 
various  theoretical  grounds  (refs.  10  and  11).  Kierstead’s  recent 
measurement  of  the  quantity  (9p/9T)^  along  the  critical  isochore  of 
He^  (ref.  12),  combined  with  Moldover^s  heat*capacity  data  (ref.  13), 
strongly  support  this  hypothesis.  In  fact,  it  appears  that  the  quantity 
(9g/9T)^  is  a  linear  function  of  T  for  P=P^  over  an  appreciable  interval 
extending  both  above  and  below  T^.  Such  a  hypothesis  also  has  a  certain 
appeal  because  it  gives  rise  to  simpler  and  more  S3mimetric  expressions 
for  a  number  of  related  properties,  and  we  will  presume  its  validity  in 
the  following  section.  It  then  becomes  natural  to  adopt  the  free  energy 
per  unit  volume,  a(p,T),  from  which  the  chemical  potential,  pressure, 
entropy  per  unit  volume,  and  heat  capacity  per  unit  volume  are  obtained 
as  follows: 


g 

=  (9a/9p) 

(19) 

p 

=  pg  -  a 

(20) 

ps 

=  -(3a/3T) 

(21) 

pc 

^  V 

=  -T(32a/9T2). 

(22) 

VI.  THE  LOCAL  APPROXIMATION  FOR  THE  ANOMALY 

At  this  stage,  it  is  instructive  to  consider  an  expression  which 
can  be  valid  only  in  the  immediate  neighborhood  of  the  critical  point — for 
reasons  that  will  be  made  clear  later.  Nevertheless,  a  great  deal  may  be 
inferred  from  the  behavior  of  certain  quantities  derived  from  it  in  this 
restricted  region.  To  simplify  the  notation,  it  is  convenient  to  introduce 
as  independent  variables 

cf  =  ?  -  1  [reduced  density  difference]  (23) 

T  =  6-1  [reduced  temperature  difference]  (24) 
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We  then  write 

for  the 

singular  part 

of  the  free  energy 

a*  = 

-RTp)|;*(a,T) 

(25) 

where 

0T  . 

2  24J 

(26) 

0)  = 

T 

52  +  T 

£  1  if  T  > 

(27) 

a  = 

a  +  6t 

(28) 

and  a,  3,  Yl  ^^^d  6  denote  adjustable  parameters.  The  magnitude  or 
’'strength”  of  the  anomaly  is  measured  by  a  and  its  range  of  influence, 
by  3;  while  Yl  controls  the  curvature  of  the  coexistence  curve  at  the 
critical  point  and  6,  the  slope  of  its  diameter  (defined  for  subcritical 
temperatures  as  the  locus  of  the  midpoints  of  its  vapor  and  liquid  branches) . 

f 

Not  only  does  formula  (26)  possess  the  desired  type  of  singularity 
at  the  critical  point  (a=0,T=0),  but  the  argument  of  the  logarithm  also 
tends  to  infinity  everywhere  along  the  locus  t  =  ~yi5^5  thus  generating 
a  continuum  of  singular  points  which  curves  downward  from  the  critical 
point,  as  depicted  in  Figure  3.  Such  a  locus  must  lie  within  the  region 
of  the  density-temperature  plane  corresponding  to  unstable  states, 
so  that  (except  for  the  critical  point)  it  is  beyond  the  reach  of  direct 
experimental  investigation.  Its  existence  is  implied,  however,  by  the 
observation  that,  along  isochores  on  either  side  of  the  critical  one, 
the  value  of  c^  rises  at  an  ever  increasing  rate  with  falling  temperature 


It  is  shown  in  reference  2  that  formula  (26)  gives  rise  to  a  logarithmic 
specific-heat  anomaly  while  at  the  same  time  the  quantities  s*,  p*, 
9p*/9p,  9p*/8T,  9^p*/9p^,  9^p*/9p9T  and  9^p5*f/9p^  all  exist  and  remain 
bounded  as  the  critical  point  is  approached  from  the  upper  half-plane 
(t>0) •  This  is  necessary  to  ensure  that  the  formulation  represents 
a  thermodynamically  stable  situation  along  the  isotherms  just  above 
the  critical  point. 
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Figure  4  ~  Local  Behavior  of  Contours  of  o) 


until  the  coexistence  curve  is  intersected  (ref.  14).  Indeed,  the 
shape  of  this  curve  may  be  largely  determined  by  the  proximity  of 
such  a  line  of  singular  points. 

Additional  insight  is  gained  by  examining  the  family  of  trajectories 
generated  when  the  quantity  w  is  held  constant.  In  particular,  o)  takes 
on  the  value  unity  along  the  diameter  of  the  coexistence  curve,  zero  on 
the  critical  isotherm,  and  negative  infinity  along  the  line  of  singular  ' 
points,  as  illustrated  in  Figure  4.  It  can  further  be  shown  that  the 
coexistence  curve  coincides  (at  least  in  the  critical  region)  with  one 
of  the  contours  of  co.  The  regular  part  of  the  free  energy  is  expanded 
in  a  Taylor  series  about  the  critical  point,  and  to  this  is  added  the 
singular  part  given  by  equation  (25).  Then,  using  equation  (19)  to  obtain 
g(a,T)  and  invoking  the  condition  for  equilibrium  of  the  two  phases, 
g  (t)  =  g  (t)  ,  we  can  find  an  expression  for  the  value  03=0)  appropriate 

V  JLi  S 

to  this  curve.  Moreover,  as  pointed  out  in  reference  2,  certain  deriv- 
atives  of  ip  (o,r)  reduce  to  polynomials  in  00,  so  that  their  values  are 
invariant  along  a  trajectory  for  which  (jj  is  constant.  Lastly,  we  note 
that  the  inclusion  of  the  term  6t  in  the  definition  (28)  for  o  allows  the 
approximation  (27)  for  o)  near  the  critical  point  to  be  carried  as  far  as 
third-order  terms  in  a  (since  t  is  of  the  same  order  as  there).  As  a 
result,  the  parabolic  configuration  of  the  trajectories  is  skewed  so  as 
to  be  S3niimetrical  about  the  rectilinear  diameter  rather  than  about  the 
critical  isochore.  [See  reference  15  for  a  somewhat  similar  treatment 
of  the  critical  region.] 

VII.  REQUIREMENTS  FOR  GLOBAL  VALIDITY 

In  actuality,  the  coexistence  curve  can  be  characterized  as 
parabolic  only  over  an  interval  extending  to  either  side  of  the  critical 
density  by  roughly  ten  percent.  The  curvature  is  so  slight  for  most 
substances,  however,  that  the  corresponding  drop  in  temperature  amounts 
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to  only  about  a  tenth  of  a  degree!  At  greater  density  deviations, 

higher-order  terms  have  a  significant  effect,  thus  obscuring  the 
i 

parabolic  region  whose  existence  has  long  been  a  matter  of  specu¬ 
lation  (ref.  17). 

To  rectify  the  situation,  equations  (26),  (27)  and  (28)  for  i|^*, 

0)  and  a,  respectively,  must  be  so  modified  that  they  apply  throughout 
the  density-temperature  plane  while  their  local  properties  remain 
intact.  In  particular,  the  contours  of  oj  must  conform  to  a  set  of 
constraints  that  ensures  proper  asymptotic  behavior  for  extreme  values 
of  density  and  temperature.  The  desired  effect  is  achieved  by  requiring 
that  the  numerator  of  the  expression  for  co  vanish  not  only  at  the  critical 
temperature  but  also  as  the  temperature  approaches  zero.  It  must  also 
vanish  when  the  density  goes  to  zero  or  becomes  indefinitely  great.  At 
the  same  time,  the  denominator  must  be  given  additional  flexibility  to 
permit  essential  nonlocal  conditions  to  be  imposed  on  the  shape  and  location 
of  the  line  of  singular  points.  However,  any  modification  has  to  be  done 
in  such  a  way  that  certain  key  features  are  retained:  the  denominator 
function  must  remain  positive  for  all  physically  realizable  states ;  the 
line  of  singular  points  (along  which  the  denominator  vanishes)  must  be 
concave  downward  with  its  maximum  ordinate  occurring  at  the  critical 
point;  and  the  conditions  w  =  1  and  9a)/9C  =  0  must  hold  along  its 
diameter,  as  before. 

These  considerations  suggest  the  following  generalization  of 
equation  (27) 


0)  =  q^0^(0^-l)/z 


(29) 


This  phenomenon  can  best  be  seen  in  Weber’s  measurements  of  the  density 
of  oxygen  in  the  critical  region  (ref.  16).  His  data  include  six  pairs 
of  saturation  densities  that  lie  within  ten  percent  of  p  ,  and  are  the 
most  accurate  available  at  the  time  of  this  writing.  ^ 
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where  z 


Ir^  +  Jr  +  K  +  (Lr  +  M)(03-l)  +  N(e 3-1)2 


(30) 


g  +  1 

Ca^  +  a  +  1 


1  when  5=0 
0  when  5  =  -l,oo 


r 


Da^  +  1 


0  when  5=0 
1  when  5  =  -1,1 


5  =  1 


0  when  ^  =  1/u 
-1  when  5=0 


(31) 

(32) 

(33) 


and  u  =  u(0)  is  a  function  of  temperature  (to  be  given  later)  whidh 
satisfies  the  condition  u(l)  -  1,  so  that  5(5,1)  =  a.  The  quantities 
C,  D,  I,  J,  K,  L,  M  and  N  are  adjustable  coefficients  subject  to  certain 
constraints.  We  assume  that  both  C  and  D  have  nonnegative  values  so  that 
the  functions  q(5)  and  r(5)  will  be  free  of  poles  on  the  real  axis  for 
all  5  >  -1.  This  avoids  the  introduction  of  extraneous  poles  and  zeroes 
into  the  definition  of  w  that  might  adversely  affect  its  properties. 

The  denominator  z  has  the  form  of  a  general  quadratic  expression 
in  the  variables  r  and  6^-1,  both  of  which  vanish  at  the  critical  point. 
[Note:  The  quantity  0^  is  chosen  here  as  one  of  the  independent  variables 
to  ensure  conformity  with  Debye’s  law  as  0-^0,  just  as  was  done  in 
Section  IV  with  the  formulation  of  the  nonsingular  part  of  the  free  energy.] 
To  show  that  the  linear  terms  of  this  expression  suffice  to  characterize 
the  region  lying  close  to  the  critical  point  (i.e.,  to  first-order  terms 
in  the  local  variables  5^  and  t) ,  we  substitute  the  following  approximations 
into  equations  (29)  and  (30) : 


q  1;  r  -  (I>fl)5^  \a\  «  1 

03  ^  2_.  =  (02+0+1)  (e«l)  3x  |t|  <<  1 


(34) 


Then,  setting  J  =  3yi/(D+1),  K  =  0  and  M  =  1,  and  neglecting  the  quad¬ 
ratic  terms,  we  find  that  equation  (29)  reduces  to  equation  (27). 
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The  variable  q  exerts  a  major  Influence  at  extreme  values  of 
density:  it  approximates  5^/C  as  ?  0  and  l/C^u  as  5  Con¬ 

sequently,  the  numerator  of  the  expression  for  w  will  vary  quad- 
ratically  with  5  at  low  densities  and  inversely  with  at  very 
high  densities,  provided  that  C  ^  0.  The  maximum  value  of  q(a) 
occurs  at  a=0,  although  the  shape  is  skewed  with  respect  to  this 
abscissa  due  to  the  presence  of  an  odd  power  of  a  in  equation  (31). 

On  the  other  hand,  the  shape  of  r(a),  which  has  a  minimum  at  a=0, 
is  S3nnmetrical  because  equation  (32)  involves  only  even  powers  of  a. 

It  follows  that  z  must  likewise  be  an  even  function  of  a. 

VIII.  THE  FORM  OF  THE  SINGULAR  LOCUS 

As  previously  noted,  the  vanishing  of  the  denominator  of  the  ex¬ 
pression  for  0)  generates  a  line  of  singular  points.  Although  the  exis¬ 
tence  of  such  a  locus  has  long  been  suspected  (ref.  18,  §  79),  its  true 
nature  remains  unknown.  We  propose  here  to  assign  plausible  attributes 
to  it  since  its  specification  is  crucial  to  the  development  of  a  globally 
valid  approximation  for  the  anomaly.  To  that  end,  the  denominator 
function  will  be  subjected  to  certain  conditions  and  changes  of  arguments. 

In  the  r,0^-plane  the  locus  z=0  is  a  conic  section,  the  form  of 
which  depends  upon  the  values  of  the  coefficients  of  equation  (30).  Its 
contour  is  significantly  altered  by  transformation  to  the  a,6^-plane,  and 
finally  to  the  5,0-plane,  to  produce  the  desired  result.  The  rational 
expressions  (31)  and  (32)  have  been  formulated  so  as  to  make  the  first 
derivatives  of  both  q  and  r  equal  to  zero  when  a  =  0.  It  is  apparent 
from  equation  (33)  that  this  condition  occurs  whenever  5^  =  1»  that  is, 
along  the  path  given  by 

5  =  5^(0)  =  l/u(0)  (35) 

which  corresponds  to  the  diameter  or  medial  axis  of  the  curve  z(a,0)=O. 

The  only  density-dependent  variables  that  enter  into  equation  (29)  are 
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q  and  r.  Hence,  the  quantity  Bu/SC  =  u(8a)/3a)  will  vanish  everywhere 
on  the  medial  axis,  as  required.  In  order  that  the  value  of  to  should 
be  unity  along  the  same  path,  we  must  have 

•  z(O,0)  =  03(03-1)  =  03-1  +  (03-1)2  (36) 

Comparing  this  with  equation  (30)  when  r  is  set  equal  to  0,  we  see  that 

K  =  0 

M  =  1  (37) 

N  =  1 

The  values  of  the  other  three  coefficients  remain  to  be  determined. 

Two  of  the  degrees  of  freedom  will  be  used  to  satisfy  constraints  on  the 
behavior  of  the  singular  locus  as  it  approaches  the  origin  in  the  5,0-plane. 
These  constraints  derive  from  the  observation  that  the  saturation  density 
as  well  as  the  vapor  pressure  of  pure  substances  decreases  fairly  rapidly 
between  the  critical  point  and  the  triple  point.  If  the  density  is  ex¬ 
pressed  as  a  power  of  the  temperature,  the  exponent  typically  varies  from 
around  6  in  the  upper  part  of  the  interval  to  considerably  more  than  12 
near  the  triple  point.  Thus,  the  vapor  branch  of  the  coexistence  curve 
approaches  the  temperature  axis  very  closely  at  low  temperatures. 

It  is  reasonable  to  suppose  that  the  line  of  singular  points  follows 
a  somewhat  similar  path,  descending  from  either  side  of  the  critical  point 
(where  it  has  a  common  tangent  with  the  coexistence  curve)  to  intercept 
the  density  axis  at  5=0  and  at  a  second  point  well  beyond  the  critical 
value,  as  shown  in  Figure  3.  Moreover,  it  will  be  assumed  that  this  locus 
is  tangent  to  the  temperature  axis  at  the  origin  and  that  its  density 
coordinate  varies  near  there  as  a  power  of  0.  We  take  the  exponent  to 
be  6  to  make  certain  that  the  resulting  path  will  always  lie  within  the 
coexistence  boundary.  These  requirements  are  satisfied  by  setting  both 
z  and  its  first  partial  derivative  with  respect  to  0^  equal  to  zero  at  the 
origin,  which  yields  the  relations 
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J 

L 


1  -  I 
1 


(38) 


The  remaining  degree  of  freedom  will  be  used  to  accommodate  to  the 
unique  behavior  of  the  substance  in  question. 

Let  us  write  the  coefficient  I  in  equation  (30)  as  0^  and  replace 
the  others  with  the  values  deduced  above.  The  denominator  function, 
arranged  according  to  powers  of  r  or  6^,  respectively,  then  becomes 


i^r^  +  (0^“0^)r  +  0^(0^-l) 

O  0 

(39a) 

^  -  (l-r)03  -  03r(l-r) 
o 

(39b) 

The  intersection  of  the  surface  generated  by  this  function  with  the  plane 
z=0  defines  the  singular  locus  in  the  r, 6 ^-coordinate  system.  The  resulting 
curve  will  be 

a)  an  ellipse  if  0^  >  1/4, 

b)  a  parabola  if  0^  =  1/4,  or 

c)  a  hyperbola  if  0^  <  1/4  (but  ^0). 

The  last  case  degenerates  into  two  intersecting  straight  lines  when  0^  =  0. 

We  shift  to  the  5,0 ^-coordinate  system  via  the  relation 

_ r _ 

^  1  +  D(l-r) 

obtained  by  inverting  equation  (32) .  This  will  always  yield  pairs  of  real 
roots  for  those  points  of  the  locus  having  r-coordinate  values  in  the  range 
0  to  1.  The  constraints  placed  on  z  predetermine  the  occurrence  of  just  one 
such  arc  which  (under  the  additional  restriction  0  >0)  transforms  into  a 


1/2 


D  >  0 


(40) 


In  the  hyperbolic  case  the  other  branch  of  the  locus  is  confined  to  the 
region  r  >  1,  Further,  if  negative  values  of  0^  are  disallowed,  it  will 
lie  entirely  within  the  fourth  quadrant  (where  the  temperature  is  negative) 
and  can  be  safely  Ignored  since  it  refers  to  a  nonexistent  state. 
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kidney-shaped  contour  with  bilateral  symmetry  about  the  medial  axis. 

It  has  a  horizontal  tangent  at  the  critical  point,  and  Intercepts  the 
density  axis  at  the  points  where  a  =  +1.  Extending  below  the  axis  are 
two  lobes  that  unite  to  form  an  arch  tangent  to  this  line  at  the  point 
where  5=0.  The  exact  shape  will,  of  course,  depend  on  the  parameters 
6^  and  D,  but  regardless  of  the  form  the  value  of  z  will  be  negative  at 
all  points  interior  to  the  curve  and  positive  elsewhere,  so  long  as  the 

aforementioned  restrictions  on  these  two  quantities  are  observed. 

\ 

Because  of  symmetry,  the  medial  axis  constitutes  a  locus  of  extreme 
values  for  the  isothermal  sections  of  the  surface  z(a,0).  Whether  they  are 
maxima  or  minima  depends  on  the  sign  of  the  second  term  in  equation  (39a), 
which  varies  with  temperature*  It  is  evident  that  the  curvature  of  the 
isotherms  undergoes  a  reversal  at  0=0^-  Above  the  transition  temperature, 
the  surface  is  characterized  by  a  single  central  trough,  whereas  below  this 
value  it  divides  into  two  troughs  separated  by  a  ridge  which  is  centered 
along  the  medial  axis.^  The  geometry  of  the  surface  in  the  critical  region 
will  be  especially  sensitive  to  this  parameter  whenever  its  value  is  near 
unity*  (Additional  analytical  results  are  given  in  Appendix  A*) 


i* 

It  is  interesting  to  compare  the  isotherms  of  z(a,0)  with  results  obtained 
by  Hemmer,  Kac  and  Uhlenbeck  (ref*  19)  for  a  one-dimensional  fluid  model  in 
which  the  effect  of  both  short-range  (repulsive)  and  long-range  (attractive) 
forces  is  taken  into  account*  Through  an  appropriate  change  of  variables, 
they  develop  a  successive-approximation  method  for  solving  the  Kac  integral 
equation  which  is  valid  in  the  critical  region.  To  zeroth  order,  it  yields 
the  van  der  Waals  equation  of  state*  The  next  approximation  includes  a  term 
representing  the  first  “quantum”  correction,  which  becomes  significant  only 
as  the  two-phase  region  is  approached.  The  solution  is  interpreted  by  analogy 
with  that  of  a  corresponding  Schrodinger  equation.  The  lowest  eigenfunction, 
which  is  symmetric,  has  the  form  of  a  fourth-degree  polynomial  [see  eq.  18  in 
ref*  19]  and  its  coefficients  exhibit  a  temperature  dependence  like  that  of 
equation  (39a)  above*  In  the  critical  region  there  is  a  smooth  transition 
from  a  function  with  a  single  minimum  to  one  with  a  double  minimum  (when  the 
temperature  falls  below  the  critical  value) *  It  also  turns  out  that  the 
coexistence  curve  is  appreciably  flatter  than  the  van  der  Waals  theory 
predicts,  but  it  still  remains  parabolic  at  the  critical  point. 
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IX.  THE  EXTENDED  APPROXIMATION  FOR  THE  ANOMALY 


The  ground  has  now  been  laid  for  the  proposed  generalization  of  i|;*, 
namely. 


In  -JL. 

4 


1 

2^ 


(41) 


in  which 


X  =  33a)/(0^-l)  =  3Bq^0^/z 


(42) 


X  =  x(?  ,6)  =  3B/(e^-l) 

m  m 


(43) 


We  also  furnish  an  explicit  formula  for  the  function  u(0)  which  fixes  the 
location  of  the  medial  axis  via  equation  (35), 


_ +  1.4 _ 

(1-0.46) (06-1)  +  2.4 


6  <  2.5 


(44) 


The  value  of  x  given  by  equation  (42)  is  such  that  it  approaches  zero  under 
all  extreme  conditions  of  density  and  temperature,  whereas  it  becomes  in¬ 
definitely  great  along  the  line  of  singular  points.  It  follows  that  the 

"X 

quantity  e  not  only  vanishes  at  the  critical  point  but  remains  negligible 
in  the  immediate  vicinity  as  long  as  x  is  sufficiently  large.  Indeed,  this 
property  may  be  taken  as  demarcating  the  approximate  extent  of  the  critical 
region.  Under  these  circumstances,  only  the  last  expression  within  the  braces 
of  equation  (41)  is  of  any  consequence,  and  it  reduces  to  the  form  of  equation 
(26),  as  required. 


Outside  the  critical  realm,  both  terms  of  equation  (41)  are  significant, 
although  their  respective  contributions  to  the  specific-heat  anomaly  vanish 

—X 

as  X  0.  By  expanding  the  quantity  ln[x/(l-e  )  ]  in  terms  of  a  power  series 
in  X, 


r2  4.  nv3  _L.  _ 1 


^  -X  2^  “  2880^ 

1-e 


(45) 
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we  can  show  that  the  first  terra  of  equation  (41)  varies  as  when  x  0, 

while  the  second  one  varies  in  accordance  with  the  behavior  of  x/x  (or  o)) 

in 

under  the  same  conditions.  Consequently,  the  contribution  of  the  first 

terra  becomes  negligible  under  all  extreme  conditions.  The  second  term  is 

also  found  to  be  negligible  except  when  6  becomes  indefinitely  great.  In 

this  particular  case,  its  value  approaches  a  negative  constant  [or  0  if 

5  =  5  (oo)  ]  with  the  result  that  a*  will  ultimately  Increase  linearly  with 
m 

temperature,  thus  contributing  to  the  entropy  but  not  to  the  specific  heat. 
Such  behavior  is  not  only  admissible  but  is,  in  fact,  necessary  in  order  to 
reproduce  the  observed  character  of  the’ anomaly. 

“X  “X 

The  quantities  e  and  1-e  that  multiply  the  first  and  second  terms, 

respectively,  of  equation  (41)  may  be  regarded  as  weighting  factors  which 

effect  a  transition  in  from  the  behavior  peculiar  to  the  critical  region 

’I' 

to  a  form  consistent  with  that  observed  elsewhere  .  A  special  situation 
exists  along  the  path  5=5^(0) •  the  two  expressions  enclosed  in  brackets 
become  identical  and  the  form  of  is  independent  of  the  weights.  Further¬ 
more,  along  supercritical  isotherms,  the  quantity  x  (and  therefore  5^*) 
attains  a  maximum  while  a*  passes  through  a  minimum  there.  As  a  result,  the 
singular  part  of  the  expression  for  the  chemical  potential  [g*  =  3a*/ 9p] 
vanishes  on  the  medial  axis . 

This  locus  will  not,  in  general,  coincide  with  the  diameter  of  the 

coexistence  curve,  although — just  below  the  critical  point — it  does  lie  very 

close  to  the  midpoints  of  the  vapor  and  liquid  branches,  and  ^  the  critical 

point  the  two  loci  do  have  a  common  tangent.  Equation  (44)  provides  a 

means  for  adjusting  the  slope  of  this  tangent,  since  6  =  “(d5^/d0)^» 

In  addition,  the  particular  form  given  to  u  has  been  chosen  to  make 

(d^5  /d0^)  =  0,  i.e.,  so  that  the  point  of  inflection  will  occur  at  0=1. 

m  c 

i" 

The  introduction  of  such  functions  permits  a  continuous  transformation 
from  an  appropriate  local  description  to  one  that  conforms  to  prescribed 
asymptotic  conditions.  This  gives  rise  to  an  essential  singularity  (in 
the  function-theoretic  sense)  at  those  points  where  the  exponent  goes  to 
infinity.  Such  a  possibility  was  originally  noted  by  Zimm  (ref.  20). 


23 


The  section  of  the  locus  inmiediately  below  the  critical  point  will  then 
approximate  a  straight  line,  in  conformity  with  the  well-known  law  of 
the  "rectilinear  diameter,"  That  the  diameter  may  be  so  characterized  is 
borne  out  by  Weber’s  precise  determination  of  the  saturation  densities  of 
oxygen  (ref.  16),  which  reveals  no  measurable  curvature  over  an  interval 
extending  several  degrees  below  T^. 

X.  CONCLUSION 

The  equation  of  state  proposed  here  is  intended  to  span  the  range 
of  states  from  gases  to  very  dense  liquids  and  from  low  temperatures  to 
moderately  high  ones.  It  has  been  formulated  to  exhibit  correct  asymp¬ 
totic  behavior  (at  least  in  those  cases  for  which  theory  has  already 
provided  a  definite  answer)  as  one  or  the  other  of  the  independent 
variables  approaches  zero  or  Infinity.  This  allows  the  flexibility 
provided  by  the  adjustable  constants  to  be  devoted  to  the  attainment  of 
close  agreement  with  observation  in  the  regions  where  accurate  measurements 
are  available,  with  some  degree  of  assurance  that  the  approximation  will 
not  yield  unreasonable  values  elsewhere. 

Most  of  the  parameters  occur  linearly  (as  coefficients  of  the  various 
terms),  thus  facilitating  a  weighted,  least-squares  determination  of  their 
values.  The  few  which  are  involved  in  a  nonlinear  manner  have  special 
significance.  They  have  been  introduced  to  characterize  behavior  peculiar 
to  certain  regions,  and  initial  estimates  of  their  values,  which  are  needed 
at  the  outset  of  the  fitting  procedure,  can  usually  be  inferred  from  pre¬ 
liminary  analysis  of  selected  data. 

Furthermore,  the  resulting  approximation  is  readily  differentiated 
to  obtain  expressions  for  such  quantities  as  pressure,  compressibility, 
enthalpy,  entropy,  specific  heat,  etc.  Its  derivatives  all  exist  and  are 
continuous  everywhere  except  on  the  line  of  singular  points.  ,  However, 
since  this  locus  lies  in  a  region  inaccessible  to  observation,  such 
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a  limitation  is  of  no  practical  concern  apart  from  the  critical  point, 
where  a  logarithmic  singularity  appears  in  the  isochoric  specific  heat. 

Finally,  we  list  the  chief  advantages  that  derive  from  adopting 
a  formulation  based  on  the  Helmholtz  function: 

1)  Permits  a  global  correlation  of  empirical  data  in  terms  of 
a  single  function  of  density  and  temperature 

2)  Avoids  piecewise  approximation  with  attendant  boundary 
matching  problems 

3)  Eliminates  occurrence  of  integrals  in  expressions  for 
evaluating  related  properties 

4)  Ensures  thermodynamic  consistency  of  all  derived  quantities. 


APPENDIX  A 

ANALYSIS  OF  BEHAVIOR  NEAR  THE  CRITICAL  POINT 


The  shape  of  the  section  z(d,l)  is  of  particular  interest  since  it 
largely  determines  how  both  the  singular  locus  and  the  coexistence  curve 
behave  near  the  critical  point.  Its  local  properties  can  be  analyzed  by 
first  expanding  r(a)  in  a  power  series, 

r  =  (IH-l)d^[l  -  Dd^  +  -  ...] 

and  then  substituting  the  result  into  equation  (39a)  with  0  set  equal 


(Al) 
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to  1.  Thus,  we  obtain 


where 


z(a,l)  =  +  ... 

2 

3Yi  =  i-^(0,l)  =  (1-0  3)  (D+1) 

^  ^  95-^  o 

3Yo  =  =  [e3(2I>fl)  -  D](]>fl) 

2  24  954  o 


(A2) 

(A3) 

(A4) 


We  know  on  physical  grounds  that  z  can  vanish  on  the  critical  isotherm 
only  at  the  critical  point,  where  a  singularity  appears.  Elsewhere  on  this 
isotherm  its  value  must  be  positive.  Therefore,  the  first  nonzero  coef¬ 
ficient  of  the  expansion  (A2)  must  have  a  positive  sign,  from  which  we  get 
the  condition  >  0.  In  addition,  we  surmise  that  y^  >  0  from  the  general 
observation  that  the  curvature  of  the  coexistence  boundary  increases  rather 
sharply  in  either  direction  from  the  critical  point.  These  two  conditions, 
taken  in  conjunction  with  equations  (A3)  and  (A4) ,  respectively,  imply  upper 
and  lower  bounds  on  the  permissible  range  of  0^,  viz.. 


03  i  1  if  Yi  ^  0 

0  3  D/(2]>fl)  If  Y  >  0 

O  2  ^ 


(A5) 


The  simultaneous  solution  of  equations  (A3)  and  (A4)  can  provide  ex¬ 
plicit  formulas  for  0^  and  D  in  terms  of  y^  and  y^.  These  may  be  written 


where 


03  =  X/(X  +  3Yy 
D  =  X  +  3Yj  -  1 

X  =  [3y2  -  3y^a  -  =  03(I>fl) 


(A6) 

(A7) 

(A8) 


In  order  that  X  should  be  real,  the  radicand  in  equation  (A8)  must  be  non¬ 
negative.  Hence,  we  have  the  inequality 

Y2  ^  “  3y^)  (A9) 
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Moreover,  A  is  taken  to  be  the  positive  root  because  only  this  case  can 
yield  a  value  for  D  greater  than  zero.  The  sum  of  the  coefficients 
and  Y2  is  also  subject  to  an  inequality.  By  adding  equations  (A3)  and 
(A4) ,  we  find 


3(Yi  +  Y2)  =  [1  +  (2e3-l)D](rH-l) 

=  1  +  [2e3(rH-l)  -  D]D  ^  1  (AlO) 

provided  that  9^  >  D/2(I>fl)  and  D  >  0  (which  do  not  contradict  any  previous 
assumptions).  Thus,  we  reach  the  conclusion  that  and  y^  cannot  both  be 
small  at  the  same  time. 

A  power  series  in  a  may  also  be  obtained  for  the  line  of  singular  points 
First,  we  set  z  equal  to  zero  in  equation  (39b)  and  solve  the  resultant  quad 
ratic  equation  for  0^  in  terms  of  r.  This  gives 

03  =  ■|•(l-r)  ±  |■[(l-r)2  +  4e3r(l-r)]^'^^  (All) 

of  which  only  the  positive  root  can  describe  a  locus  in  the  first  quadrant 
of  the  r,0^“-plane.  Then,  after  substituting  expression  (Al)  for  r(a),  we 
form  the  series  expansion  for  0^  together  with  its  cube  root.  The  final 
results  are 


16^(5)  ^  ~  3Yj5^  -  3(y2  -  AYi)a'+  +  ...  (A12) 

^  -  AYi  +  Yj)a'+  +  ...  (A13) 

As  expected,  the  coefficient  of  the  quadratic  term,  which  measures  the 
downward  curvature  of  the  locus  at  the  critical  point,  agrees  with  the 
local  approximation  considered  in  Section  VI.  It  is  also  evident  that, 
if  Yi  is  small,  the  quartic  term  cannot  have  a  sizable  coefficient  unless 
Y2  Is  fairly  large. 
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